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ABSTRACT 

CN , 

Context. We present a detailed study of localised magnetohydrodynamical (MHD) instabilities occuring in two- 
dimensional magnetized accretion disks. 

Aims. We model axisymmetric MHD disk tori, and solve the equations governing a two-dimensional magnetized accre- 
tion disk equilibrium and linear wave modes about this equilibrium. We show the existence of novel MHD instabilities 
in these two-dimensional equilibria which do not occur in an accretion disk in the cylindrical limit. 
Methods. The disk equilibria are numerically computed by the FINESSE code. The stability of accretion disks is in- 
vestigated analytically as well as numerically. We use the PHOENIX code to compute all the waves and instabilities 
accessible to the computed disk equilibrium. 

Results. We concentrate on strongly magnetized disks and sub-Keplerian rotation in a large part of the disk. These 
disk equilibria show that the thermal pressure of the disk can only decrease outwards if there is a strong gravitational 
potential. Our theoretical stability analysis shows that convective continuum instabilities can only appear if the density 
contours coincide with the poloidal magnetic flux contours. Our numerical results confirm and complement this theo- 
retical analysis. Furthermore, these results show that the influence of gravity can either be stabilizing or destabilizing 
on this new kind of MHD instability. In the likely case of a non-constant density, the height of the disk should exceed 
S—j \ a threshold before this type of instability can play a role. 

Conclusions. This localised MHD instability provides an ideal, linear route to MHD turbulence in strongly magnetized 
accretion disk tori. 

Key words. Accretion, accretion disks - Instabilities - Magnetohydrodynamics (MHD) - Plasmas 

03 ' 1- Introduction 

_ _ i 

Accretion disks are very common objects in astrophysics. These objects can be found in stellar systems as well as in 
active galactic nuclei (AGN). In the first case, the disk accretes matter onto a protostar (young stellar object (YSO)), 
white dwarf, neutron star, or a black hole. The typical size of these disks is a few hundred astronomical units (AU). A 
massive black hole is the central object in the case of an accretion disk in AGN. Here, the typical size of the disk is a 
hundred parsecs (pc). 

Much research on accretion disk dynamics focuses on occuring accretion processes in an MHD framework. This is 
done linearly as well as non-linearly. The linear studies aim to understand the drivers of the accretion mechanism. This 
can be done by looking at waves and instabilities about a disk equilibrium in the cylindrical limit (see e.g. Keppens et 
al. 120021 Blokland et al. 120051 van der Swaluw et al. 2003 and references therein). In this limit the disk equilibrium is 
essentially one-dimensional and one can use a self-similar model like the one of Spruit et al. (|1987p . When instabilities 
are found, one must compute the resulting evolution non-linearly to see if these instabilities give rise to turbulence. 
This turbulence may be the source of an enhanced effective viscosity mechanism which causes an increased transport 
of angular momentum outwards, as needed for accretion (Shakura & Sunyaev[T973). Balbus & Hawley (|199ip discussed 
that the MHD turbulence resulting from the magneto-rotational instability (Velikhov 1959] and Chandrasekhar 1960) 
could provide this angular momentum transport. 

In the last years, global two or even three-dimensional magnetohydrodynamical (MHD) simulations of accretion 
disks have been performed (see for example Matsumoto & Shibata [19971 and Armitage I1998[) . In these simulations, the 
non-linear dynamics is usually interpreted as a direct consequence of the magneto-rotational instability (see for ex- 
ample Hawley et al. 2001) leading to angular momentum transport outwards. Another possibility for the initial and 
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evolving dynamics is that both the convective and the magneto-rotational instability play an important role in the 
angular momentum transport (see for example Igumenshchev et al. 2003). In the latter case, these disks are known as 
convection-dominated accretion flows (CDAF). However, especially those global simulations starting from initial axisym- 
metric accretion tori (see for example Hawley 20tJD]) are at best loosely connected to the linear stability studies. A detailed 
catalogue of the MHD spectrum of such two-dimensional disks is completely lacking. Moreover, the usual identification 
of the cause of the turbulent dynamics with the linear magneto-rotational or convective instability is mostly based on 
the extrapolation of the linear stability analysis from a one-dimensional to a two or three-dimensional accretion disk. 
Such extrapolating ignores the fact that if one performs a detailed MHD stability study of a two-dimensional accretion 
disk, one may find many new types of instability that all could lead to effective angular momentum transport. A recent 
example of such new, poloidal flow-driven type of instability is presented by Goedbloed et al. (|2004ap . where the Trans- 
slow Alfven Continuum (TSAC) mode is shown to occur inside a disk with toroidal and super-slow poloidal flow in the 
presence of a strong gravitational potential. 

This paper has two aims. The first one is to present the equations and the numerical solutions of two-dimensional MHD 
disk equilibria with toroidal flows, but without poloidal flows. This case is fundamentally different from the previously 
discussed one of combined poloidal and toroidal flows (Goedbloed et al. 2004a) because the equilibria are described by 
quite different flux functions. By considering an axisymmetric equilibrium, we are able to model a thick accretion disk 
without any further approximations. All the equilibria presented below are computed using the code FINESSE (Bclicn 
et al. I2002p . The second aim is to present a detailed stability analysis of these two-dimensional equilibria. The analysis 
is done theoretically as well as numerically. For the stability computations we have used the recently published spectral 
code PHOENIX (Blokland et al. [2006). To our knowlegde, this kind of analysis for MHD disk equilibria with purely 
toroidal flow has never been presented in the astrophysical literature until now. 

We will present a sample of disk equilibria where the disk plasma typically varies from strongly magnetized up to close 
to equipartition. Also, the rotational speed of the plasma varies from sub-Keplerian up to Keplerian. The theoretical 
part of the spectral analysis of these equilibria shows that the convective continuum instabilities may appear inside the 
disk. These instabilities are part of the continuous branches that exist in the MHD spectrum of the linear eigenmodes of 
the disk tori and are localised on magnetic flux surfaces. We derive a stability criterion for this instability. This criterion 
looks similar to the Schwarzschild criterion. However, our criterion governs convective instabilities along the poloidal 
magnetic field lines while the Schwarzschild criterion applies in the direction perpendicular to the poloidal magnetic field 
lines. This theoretical analysis has been verified by our numerical stability calculations. 

The paper is organised as follows. In Sect. [23 we present the equations which govern a two-dimensional accretion disk 
equilibrium. In Sect. [31 we recall the essential elements from spectral theory of MHD waves and instabilities, such as 
the Frieman and Rotenberg formalism (Frieman & Rotenberg |1960p . straight field line coordinates and representation. 
In Sect. [H these are used to derive the equations for the continuous MHD spectrum and the stability criterion for the 
convective instability. The numerical codes FINESSE and PHOENIX are discussed in Sect.O In Sect. [SI we present our 
numerical results on the disk equilibria and their stability analysis. Finally, in Sect. we summarise and present our 
conclusions. 



2. Accretion disk equilibrium 

We consider an axisymmetric accretion disk. Because of this symmetry, cylindrical coordinates (R, Z, ip) are used and 
the equilibrium quantities only depend on the poloidal coordinates R and Z. The disk equilibrium is modeled by the 
ideal MHD equations, 

p— = -pvVv-Vj)+jxB-/)V$, (1) 
at 

dp 

— = -v • Vp - 7pV • v, (2) 
— =Vx(vxB), (3) 

f = -V.(pv), (4) 

where p, p, v, B, $ and 7 are the density, pressure, velocity, magnetic field, gravitational potential, and the ratio of the 
specific heats, respectively. The current density j = V x B and the equation V • B = has to be satisfied. Furthermore, 
the disk equilibrium is assumed to be time- independent. The relation between the density and the thermal pressure is 
expressed by the ideal gas law: p = pT. From the equations V • B = and V • j = it follows that the magnetic field and 
the current density can be written as 

1 

R~ 
1 

J = -R f 

respectively. Here, 2irip is the poloidal flux and the poloidal stream function I = RB V . We restrict ourselves to disk 
equilibria with purely toroidal flow, 

v = v tp e tp = Rfle v . (7) 



B = — e v x + B^, (5) 
xVJ + j'^, (6) 
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In this case, equations © and (U]) are trivially satisfied. The angular velocity is related to the electric field, 

E = v x B = QVip, (8) 

and making use of the induction equation ([3]) it can be easily shown that O = Q(ip). 

Three different projections can be applied on the momentum equation fT]). The projections are in the toroidal direction 
and in the poloidal plane parallel and perpendicular to the poloidal magnetic field lines. From the toroidal projection 
one can show that the poloidal stream function is a flux function, i.e. / = I(ip). The projection parallel to the poloidal 
magnetic field results in two equations, 



dp 
OR 

dp 
dZ 



tp— const 



ip— const 



<9$\ 



p m 2 - — 



dRJ ' 



<9$ 



(9) 



which have to be satisfied simultaneously. From these two equations we conclude that the pressure p = p{ip] Ri Z). The 
last projection, perpendicular to the poloidal magnetic field, leads to the extended Grad-Shafranov equation, 



i? 2 V • 



1 

R 2 



Wip 



dJ 

d^ 



-I- 



R 



2 dp 
dip' 



(10) 



The two equations parallel to the poloidal magnetic field ([9]) can be solved analytically under the extra assumption 
that cither the temperature T, or the density p, or the entropy S = pp^ 1 is a flux function. The assumption that the 
temperature is a flux function can be justified due to the high thermal conductivity along the magnetic field lines. This 
is true at least up to the transport time scale, which is long compared to the Alfven time. The resulting pressure can 
then be written as 

,2 D 2x A , M *(R,Z)' 



p(ip;R,Z) =_p o (-0)exp 



{if - Rt)Mi>) 



T(V>) 



(11) 



where At = fl 2 /(2T) is a flux function, po is the pressure for a static pure Grad-Shafranov equilibrium without gravity, 
and R is the geometric axis of the accretion disk. The extended Grad-Shafranov equation (fTU|) reduces to 



R 2 V 



R 2 



dJ 

'dtp 



-I- 



R 2 



dpp 

dip 

x exp 



Po 



/ R 2 R 2N dA T 



$ dT 

T 2 Ihp 



, R — Rn) At — 



T 



(12) 



Another possibility, on MHD time scales, is to assume that the density is a flux function. In this case, the pressure 
reads 

1 + (R 2 -R 2 )A P W 



p(ip;R,Z) =p (tp) 



(13) 



where the quasi-temperature T p = po/p and A p = fl 2 /(2T p ). Under this assumption, the extended Grad-Shafranov 
equation (|10p can be written as 



R 2 V 



R 2 



dpp 

dip 



- Po 



(R 2 R 2^A P $ dT p 



1 + (R 2 - Rl)K p 



fa 



1 + (R 2 - R 2 )K P - — 



(14) 



The final option is to assume that the entropy S is a flux function. The advantage of this assumption is that it allows 
for a natural extension to include both toroidal and poloidol flows in the equilibrium, where the entropy has to be a flux 
function (Zehrfeld et al. 119721 Hameiri[1983). In this case, the pressure reads 



p(iP;R,Z)=p (iP){l + 



7-1 
7 



(R 2 - R 2 )A s (iP) 



<HR, z) 



and the extended Grad-Shafranov equation (flTJ)l reduces to 



7/(7-1) 



(15) 



R 2 V 



R 2 



Vip 



dpp 
dip " 

1 + 



-Po 

7-1 

7 



(R 2 - R, 



dAc 



4> dT s 



01 dip T 2 dip 



7-1 



(j^-iZg)As-l. 

is 



7/(7-1) 



(R 2 



R 2 o)^s 



(16) 
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where the quasi-temperature Tg = ffpp - 1 and A5 = £l 2 /(2Ts). The extended Grad-Shafranov equation has been used 
previously by van der Hoist et al. (|2000b[) for the first two cases without gravity. 

The density can be easily derived for all three assumptions by inserting the corresponding equation for the pressure 
into the momentum equations parallel to the poloidal magnetic field lines ([9]) . The resulting density is 



(R 2 - R 2 ) A T - - 



p(iP;R,Z) = poO) x < 



1 



7-1 

7 



R 2 



Rl) As - ± 



1/(7-1) 



for T = T{ip) 
for p — p(ip) 

for S = S(i/>) 



(17) 



The flux function po corresponds to the density of a static equilibrium without gravity. All these cases with the inclusion 
of gravity have been discussed in the appendix to Blokland et al. (20011). 



3. Spectral formulation 

3.1. Frieman-Rotenberg formalism 

For the investigation of stability properties of stationary MHD equilibria, the formalism by Frieman and Rotenberg (1960) 
has been used. They derived from the linearised MHD equations, one second order differential equation for the Lagrangian 
displacement vector field 

F(O-2pv.V§-4f=0, (18) 

where the force operator F(£) is 

F(0 = F static (0 + V$V • (p£) + V • [p£v • Vv - pvv • V£] . (19) 

Here, 

F sta tic(£) = -Vn + B- VQ + Q VB (20) 

is the force operator for static equilibria without gravity, which has been derived by Bernstein et al. (|1958[) . Here, the 
Eulcrian perturbation of the total pressure, 

n = - 7P V-£-£ • Vp + B-Q, (21) 
and the Eulerian perturbation of the magnetic field, 

Q = V x (£ x B) . (22) 

The time-dependence for the displacement field is assumed to be an exponential one with normal mode frequencies u>, 
£ = £ exp(— Vjjt). Using this assumption, the Frieman-Rotenberg equation can be written as 

F(£) + 2ipwv • V£ + plo 2 ^ = 0. (23) 

In the derivations below, the toroidal symmetry of the accretion disk equilibrium is exploited. This is done by writing 
the toroidal dependence of the displacement field as £ ~ exp(irup), where n is the toroidal mode number. 

3.2. Straight field line coordinates 

The analytical and numerical analysis of MHD waves and instabilities are done in 'straight field line' coordinates. To 
convert from cylindrical (R, Z, <p) to straight field line coordinates (x 1 = ip, x 2 = i9, x 3 = ip) one needs the metric tensor 
and the Jacobian associated with the non-orthogonal coordinates in which the equilibrium field lines appear straight. 
This is standard practice in MHD spectral studies for laboratory tokamak plasmas. The metric elements gij and the 
Jacobian J are 

dr dr 



V.<" • V.r'. 9ij = — .—, J=(VV>X W-V(p) \ (24) 

respectively. Here, the poloidal angle ■d is constructed such that the magnetic field lines are straight in the (1?, <p)-plane. 
The slope of these lines is a flux function: 



dip 
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where q is the safety factor. Furthermore, in straight field coordinates the poloidal and toroidal curvature of the magnetic 
surfaces are 

d^) JB § , (26) 

922 J 

—d#) R, (27) 

522 / 

respectively. Here, n = Vip/\Vip\ and t = H.g/B.g, where B$ is the poloidal magnetic field. It is important to realise that 
the straight field coordinates can only be constructed when the solution ?/>(-R, Z) has been computed from the extended 
Grad-Shafranov equation (fTU]) . 



-n ■ (t • Vt) 



R 



-n • (e v ■ Ve v ) = 13$ d, 



3.3. Field line projection and representation 

On each flux surface the distinction between two wave directions can be made: one parallel and the other one perpendicular 
to the magnetic field. These directions correspond to the short wavelength limit of slow and Alfven waves in cylindrical 
geometry. Hence, it is useful to exploit a projection based on the magnetic surfaces and field lines using the triad of unit 
vectors, 

Vib B 
Using these unit vectors, the components of the displacement field £ can be written as 



X = RB#£ • n, Y = i- 
and the projected gradient operators become 

T> = -J—n • V 



B 
RB, 



RBtf 

Q = —iRBflBn ■ V 
T = -iB • V 



a #12 o 
= O-d, Of), 

922 



-iJ 
1 



- nB#, 



-i nq 

= y a * + y 



(29) 

(30) 
(31) 
(32) 



The second equality in the expressions for T and Q only holds when these operators act on a component of the displace- 
ment field 

Using the straight field line coordinates and the projections (l29l) the Frieman-Rotenberg equation (JTSj) can be written 

as 

fX\ 



(A + R)x + 2pwfiCx - puj 2 Bx = 0, 




(33) 



where u>{tp) = oj — nfl(ip) is the Doppler shifted eigenfrequency. Here, A and B are 3x3 matrix operators which are 
also present in the case of a static equilibrium without gravity. The matrix operators C and R contain the elements due 
to the toroidal rotation of the plasma and of an external gravitational potential. The matrix elements of A and B are 



An 



A 22 

A 23 
A 3 i 

A 32 
A33 



V 



R 



R 



1 B v K t 
J B,i 



1 



-v 1P g-^ -vg + 2 



A12 

A13 = -V'ypJ 7 , 

a 2 i = ^G^vU + g-v ] j + 2 



TlBtfKp B^Kt i 

R B$ J 

nBtfK 



do 



B 2 J 

j^QipQ^ 
1 



J 



G £2 G 



R 



.R 2 B 2 



T- 



B 2 



B 2 



: FjpJ 7 , 



J B$ 



(34) 
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and 



Here, the operators 



Bn 



B?Bl 



B 



22 



912 \ 
922 ) 



R 2 B 2 
B 2 ' 



B 



33 



do 



B . 

912 
922 



(35) 



(36) 



are related to the normal gradient operator T>, which has been introduced by Goedbloed (|1997p for the spectral analysis 
of static tokamak equilibria. The square brackets in the expressions (|34| for the matrix elements of A and the gradient 
operators defined by (f3T>|) indicate that the differential operator only acts on the term inside the bracket. This notation is 
also used in the expressions below. The matrices A and B were derived by Goedbloed ( |1975|, I1997[) for static equilibria. 
The matrices R and C enter if there is toroidal flow and/or external gravity. The expression for these matrices are 



Jt£> 

_P_ 

B 2 



Rn t 
Bo 



vn 2 



UpjVU 



gx 



(iDIp - 

B A 



(iVp + XT) p 



\ 



TX 



I 

B 2 



J 

d$p 



d$p 



I 

B 2 



J 



m/x 



(37) 





Rn t B 2 



Rn t B 2 R Kt I \ 



B$ 
Rn t 



B 2 
I 



Bv B 2 

R 



\ B# R 2 



where 



I? 



A 



J 

RK t 

B§ 



dsR 



9 2 



J 



J 



tOaR 



B# R 2 
jd#R 





(38) 



P J 



■ Vp, 



Vt 2 - £>$. 



(39) 
(40) 



Notice that the term p, represents the pressure variation on a flux surface. The matrix C represents the Coriolis effect 
due to the rotation while the matrix R contains rotational as well as gravitational effects. The case without an external 
gravitational potential has been discussed by van der Hoist et al. (|2000b[) . 

Two kinds of cylindrical limits can be obtained from the spectral equation (|33p . The first one is the infinitely slender 
torus, meaning that the radial position of accretion disk is taken to be at infinity. In that case the matrices R and 
C disappear. This is due to the fact that all equilibrium quantities will be become independent of the angle the 
gravitational potential at infinity is zero, and the toroidal curvature Kt becomes zero. The flow enters only as a Doppler 
shift, —nfl^), in the Doppler shifted eigenfrequency w(^). 

The other limit is the cylindrical limit of a thin (\Z\/R <§; 1) slice of plasma at the equatorial plane of the accretion 
disk. In this region all equilibrium quantities only depend on the radius R, approximately. The matrices R and C do not 
disappear as in the previous limit. Also the toroidal curvature Kt is non-zero but instead the poloidal curvature k p = 0. 
In this case the spectral equation (f3"3")l reduces to the matrix equation for a one-dimensional accretion disk presented by 
Keppens et al. PUTS)) . 



4. Continuous MHD spectrum 

In the previous section, the spectral equation (|33[) governing all MHD waves and instabilities in axisymmetric accretion 
disks has been derived. This equation is the starting point to all following MHD spectral computations. In particular, 
we can derive the equations for the continuous MHD spectrum by considering modes localised on a particular flux 
surface ip = ipQ. Van der Hoist et al. (|2000ap showed that the toroidal flow can drive the continuous spectrum unstable 
or overstable for non-gravitating plasma tori. Here, we show that the combination of toroidal flow and gravity can also 
drive the MHD continua overstable. 



4.1. General formalism 

To derive the equations for the continuous spectrum, the normal derivative (8/ d^p) of the eigenfunctions is considered 
to be large compared to the eigenfunctions themselves and the Doppler shifted eigenfrequency uj{4>) is assumed finite. In 
this case, the first row of the spectral equation (]33[) can be solved approximately, 



J 7p + B 2 \ B 2 J -fp + B 2 1P + B 2 \B 2 J V ' 
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Notice, that this solution implies that dX/dip, Y, and Z are of the same order, but more importantly that X is small 
compared to Y and Z . This means that the continuum modes are mainly tangential to a particular flux surface. 

Inserting the expression (|41[) into the second and third row of the spectral equation (|33[) results in an eigenvalue 
problem which is independent of the normal derivative. Hence, the reduced problem becomes non-singular. Exploiting 
this property, we can write the projected displacement field components Y and Z as follows: 



Y 
Z 



Sty 



(42) 



where 8{ip — ipo) is the Dirac delta function. Here, Y and Z have been split in an improper (^-dependence) and proper 
part (^-dependence). This kind of splitting has been introduced by Goedbloed (1975) for static equilibria without gravity. 
The reduced, non-singular eigenvalue problem becomes 



(a + r + 2p<jj£lc - puj 2 b) 



0. 



(43) 



where 



'^—2-^ + 4 ——(RB^k ) -2i ——j (RB#Kg) F 

B z 7p + B z 1P + B 

ypB -fpB 2 

2\J- —^77 (RBdKg) T —T 



7P + B 2 



jp + B 2 ' 



(44) 



ol 2 

4— rf ~~„„ (RB$K g ) + ^N 2 



B 2 jp + B 2 



U- 



PP B f>2 T 1 , nJN 2 



B 2 m ' po1 " 7 p + B 2 " " B 2 

Plj P2 , n T]\l2 _p>2 

\ B 2 1P +B 2 + ^ iJ Vpol D [JB 2 " U 1P + B 2 



(45) 



d#R 



jd#R 




(46) 



/ 




In these expressions k„ is the geodesic curvature, 



(47) 



JRBvB 



:d#B. 



(48) 



Here, k = b • Vb is the field line curvature and N 2 n po] is the magnetically modified Brunt- Vaisala frequency projected 
on a flux surface, 



V 2 , = JL 

m,pol g2 



1 a 

-r d vp r~52 

J p 7p + B A 



Vp- 



"B^ • Vp 






[ pB _ 




_pB ' I 



p 



jp + B 2 



Vp 



(49) 



This magnetically modified Brunt- Vaisala frequency is similar as the one presented by van der Hoist et al. (2000b) for 
tokamaks with purely toroidal flow and as the one for static gravitational plasmas published by Poedts et al. (1985) and 
Belien et al. (|1997p . Furthermore, notice that the poloidal variation of the pressure presented by the matrix r shows up 
in all its matrix elements. This destroys the diagonal dominance of the matrix a, which corresponds to the separation of 
Alfven and slow continuum modes. Similar as the matrix C in spectral equation (|33|) . the matrix c represents the Coriolis 
effect. 

The non-singular eigenvalue problem (|43[) is solved for poloidally periodic boundary conditions for the eigenfunc- 
tions r/, Solving this problem for a given toroidal mode number n on each flux surface separately results in a set of 
discrete eigenvalues. All these discrete sets together map out the continuous spectra. In the MHD formulation with the 
primitive variables p, v, p, B, an additional entropy continuum nVL is found. This Eulcrian entropy continuum does not 
couple to any of the other continua or stable and unstable modes (Goedbloed et al. I2004b|) . 



s 
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4.2. Spectral properties and stability criterion 

In this subsection a stability criterium for the continua will be derived. For this derivation we need to construct a Hilbcrt 
space with appropriate inner product. The parallel gradient operator T is a Hermitian operator under the inner product 



Jv* {Tw) dtf 



(50) 



for poloidally periodic functions v and w. Inspired by this property a Hilbert space can be defined with the inner product 



(v,w) 



v* ■ w Jd$. 



(51) 



In this Hilbert space, the matrix operators a, b, c, and r are Hermitian operators. By setting v and w equal to the 
eigenfunctions (rj, () T , one can derive a quadratic polynomial for the eigenfrequency to from the spectral equation (|4"3")l 
for the continua. 

auj 2 - 2bZ) - c = 0. (52) 

In this equation, the coefficients are 



2b = 2i i pSl 




B,B,<)K„ pu f I 



Jdi?, 



(53) 
(54) 

(55) 



where N^ v pol is the Brunt- Vaisala, frequency projected on a flux surface, 



Ni 



= JL 

BV,pol 



1 

Tp 

■ Vp 



dap 



P B 



P B 



_ EE 

ip 

B^ 
~p~B' 

B,j 



vs 



IP 



-Vp 



(56) 



1 BS 



Note that the coefficient a is always non zero but, more importantly, that the coefficients a, b, c are real due to the 
Hermitian property of the inner product. Solutions of this polynomial are u> = (b ± \fb 2 + ac)/a. This means that if 
b 2 + ac > 0, the solutions are waves with real frequencies. But if b 2 + ac < 0, the solutions contain a damped stable 
wave and an overstable mode. For the complex solution w one derives its absolute value \u\ = \f—cja by taking the 
linear combination of the polynomial with its conjugate. Furthermore, Re(w) = ntt + b/a. It can easily be shown that 
the continuum is stable if c > 0. The coefficient c is always equal or greater than zero if 



iV, 



BV,pol 



> 



(57) 



is satisfied everywhere in the plasma. From this stability criterion, we conclude that equilibria with 5* = S(ip) and equi- 
libria with T = T(tjj), 7 > 1 are always stable. Equilibria where the density is a flux function violate this criterion, which 
may result in the appearance of damped stable and overstable modes. Notice that the stability criterion (|57|) is analo- 
gous to the Schwarzschild criterion for convective instability. The only important difference is that in the Schwarzschild 
criterion one deals with normal derivatives while this criterion deals with tangential ones. This has also been noticed by 
Hcllsten & Spies (fi"979l) . Hameiri (fT983|) . and Poedts et al. (fT985]) . 



5. Numerical codes 

For the stability analysis of axisymmetric accretion tori two numerical codes, the equilibrium code FINESSE (Belien 
et al. I2002|) and the spectral code PHOENIX (Blokland et al. I2006[) . have been used. In the next two subsections their 
algorithmic details will be briefly discussed. 
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Fig. 1. An example of the geometry of an accretion disk with a circular cross-section. Here, Ro and a are the major 
radius of the geometric axis and the radius of the last closed flux surface, respectively 

5.1. The equilibrium code FINESSE 

The FINESSE code developed by Belien et al. (|2002|) can compute a stationary axisymmetric, gravitating MHD equi- 
librium described by the extended Grad-Shafranov equations (TT21) . (fl"4")) . or (flB)) for a given poloidal cross-section and 
inverse aspect ratio e = a/ Ro- Here, a and Ro are the minor radius of the last closed flux surface and the major radius 
of the geometric axis, respectively. Fig. [T] shows an example of the geometry of an accretion disk with a circular cross- 
section as seen in a poloidal plane with the central gravitational object at the origin. The extended Grad-Shafranov 
equation has been discretized using a finite element method in combination with the standard Galerkin method. The 
elements used are isoparametric bicubic Hermite ones. These bicubic elements ensure that the computed solution has the 
desired high accuracy (fourth order) needed for the stability analysis. FINESSE solves the elliptic PDE problem using 
the Picard iteration scheme. The imposed boundary conditions assume that the fixed boundary represents the last closed 
flux surface. The same FINESSE code can actually handle the more general case of stationary MHD equilibrium with 
non-vanishing poloidal flows as well, but here we restrict our discussion to equilibria with purely toroidal flows. 

5.2. The spectral code PHOENIX 

The spectral code PHOENIX developed by van der Hoist (Blokland et al. I2006P is used for the stability analysis of 
a two dimensional accretion disk. PHOENIX can compute the complete MHD spectrum for a given two dimensional 
equilibrium, toroidal mode number n, and a range of poloidal mode numbers m. A range of poloidal mode numbers is 
needed due to the mode coupling caused by allowing for a non-circular cross-section, gravitational stratification, and 
the toroidal geometry of the accretion disk. PHOENIX does not directly solve the spectral equation (|3"5)) . but instead 
employs the linearised version of the full set of MHD Eqs. ©-((I]). These linearised equations are then discretized using 
a finite element method in the normal ^-direction, a spectral method in the poloidal direction and a standard Galerkin 
method is used to obtain a linear generalised eigenvalue problem for the eigenfrequencies uj. The elements used in the ip~ 
direction for the perturbed quantities are a combination of quadratic and cubic Hermite elements to prevent the creation 
of spurious eigenvalues fRappaz I1977p . This results in a non-Hermitian eigenvalue problem, which is solved using the 
Jacobi-Davidson method (Sleijpen & van der Vorst [TOM]) . The boundary conditions used for the perturbations are those 
of a perfect conducting wall. Strictly speaking, these conditions can not be applied to waves and instabilities that appear 
inside an accretion disk. However, if the wave phenomena are sufficiently localised or not significantly affected by the 
position of the wall, these boundary conditions have hardly any influence. This is particularly true for the continuous 
branches of the MHD spectrum on which we concentrate in what follows. These modes are purely localised on a flux 
surface, and are not affected by the wall. Like FINESSE, PHOENIX can handle poloidal flows as well, but we have 
exploit it here for equilibria with toroidal flow only. 

The continuous spectrum is computed by a less expensive method introduced by Poedts & Schwarz (1993]). In this 
method, on each individual flux surface the cubic Hermite and the quadratic elements are replaced by log(e) and 1/e, 
respectively. Due to this replacement, the singular behaviour of the continuous spectrum has been approximated by 
the perturbed quantities. This singular behaviour has been described by Goedbloed (|1975[) and Pao (|1975|) for static, 
axisymmetric and toroidal plasma. Recently, the analogous treatment has been used by Goedbloed et al. (2004a) for 
plasmas with toroidal and poloidal flows and gravity. The resulting eigenvalue problem is then a small algebraic problem 
of a size set by the range in poloidal mode numbers m on each flux surface, and is then solved using a direct QR method 
and a scan over all flux surfaces. This yields detailed information on all MHD continua. 

6. Accretion disks with density as a flux function 

We present our numerical results on disk equilibria and stability using the codes described in the previous section. For the 
equilibrium computations the ratio of the specific heats 7 = 5/3 and the gravitational potential is an external Newtonian 
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one, i.e. $ = —GM*/^R 2 + Z 2 , where G and M* are the gravitational constant and the mass of the central object, 
respectively. Furthermore, we assume that the density is a flux function for all equilibrium computations. From our 
theoretical analysis above, this is the relevant case which may give rise to damped stable and overstable continuum mode 
pairs due to the presence of toroidal flow and gravity. 

6.1. Thin accretion disk with constant density 

As the first case, we take an accretion disk with constant density as also presented by Blokland et al. (|2006|) . The 
equilibrium is described by the following flux functions: 

J 2 0) =A(l- 0.0385V- + 0.02^ 2 + 0.00045^), p(ip) = 1, 

p (il>) = AB(l-0.9i>), = C(l - 0.9V- 2 ), '" M 

where the values of the coefficients A, B, and C are given by 112, 0.01, and 0.1. Furthermore, the considered poloidal 
cross-section is circular and the inverse aspect ratio e = 0.1. In these calculations, GAf* = 1 after scaling with respect to 
a reference density, the vacuum magnetic field, and the minor radius of the accretion disk. Fig. [5] shows the computed 
two-dimensional pressure and plasma beta /3(= 2p/ B 2 ) profile. The pressure decreases monotonically outwards due to 
the fact that the term proportional to the gravitational potential $ dominates in the pressure equation (|13|) . In contrast, 
the plasma beta increases monotonically outwards from 0.63 at the inner part to 0.76 at the outer part of the disk. The 
ratio i^/uKepier = [0.085, 1.049], where i>Kepier is the Keplerian velocity. 



Pressure (10"') 
3.22 3.30 3.38 3.46 3.54 3.62 3.70 3.78 3.86 3.94 




-1.0 -0.5 0.0 0.5 1.0 



Fig. 2. The two-dimensional pressure (gray-scale) and plasma beta f3 = 2p/B 2 (contours) profile for an accretion disk 
with constant density. The cross-section is circular and the inverse aspect ratio e = 0.1. The central object is to the left 
of the figure. 

The computation of the continuous spectrum of this equilibrium is done for toroidal mode number n = — 1 and 
poloidal mode numbers m = [—3,5]. Fig. [3] shows the sub-spectrum of the MHD continua. The plots clearly indicate the 
existence of stable and overstable modes in the MHD continua. The plots (a) and (b) show the real and imaginary part 
of the eigenfrequencies ui that belong to the continuous spectrum as a function of the radial flux coordinate s = yftji. 
The frequencies of the continua in the complex plane are shown in the plot (c). It should be noted that it is difficult 
to label the different branches of the MHD continua by a "dominant" poloidal mode number m. This could be done by 
making use of analytical expressions for an infinitely slender torus (e — ► 0), but then the gravity of the central object 
would be neglected. In our case, these expressions can hardly be used to identify the different branches of the MHD 
continua, because gravity plays a dominant role in the equilibrium as well as in the spectral analysis. Notice that in plot 
(a) two branches merge at s sa 0.59, while at the same location the imaginary part of the eigenfrequency starts to become 
non-zero. Intricate mode couplings thus give rise to the emergence of pairs of overstable and damped stable modes. 

Next, we investigate the influence of the external gravitational potential on the overstable modes presented in the 
previous figure. This is done by varying the mass of the central object GM*. The result is shown in Fig. [4] where 
the most overstable mode of the MHD continua is plotted as function of the gravitational potential on the magnetic 
axis i?M- Notice that each point in this plot corresponds to a different equilibrium calculated with FINESSE. PHOENIX 
is used to analyse the stability of the computed equilibrium. The gravitational potential is stabilizing from to -0.01 and 
destabilizing from -0.01 to lower values. The zero value corresponds to a toroidally rotating, non-gravitating "tokamak" 
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plasma. Its instability is in agreement with the results by van der Hoist et al. (2000a). Fig. [5]shows the projected Brunt- 
Vaisala frequency for the points (a), (b), and (c) indicated in Fig. [U Plot (a) shows that the minimum value of this 
frequency is reached inside the plasma, while plot (c) shows that this value is reached at the edge of the accretion disk. 
Point (b) is the transition from the minimum value inside to the edge of the disk as seen in plot (b) . 

6.2. Thin accretion disks with non-constant density 

To model a more realistic accretion disk equilibrium, we allow a density gradient across the poloidal magnetic surfaces. 
All other flux functions are kept as in the previous case, i.e. 

I 2 (ip) = A(l - 0.0385?A + 0.02?/; 2 + 0.00045</> 3 ), p(ip) = 1 - 0.4-0, 

p (ip) = AB(1 - 0.9V>), Oty>) = C(l - O.90 2 ), ^ 

where the coefficients A, B, C are given by 91.5, 0.01, and 0.1. Also in this case, the poloidal cross-section is circular, the 
inverse aspect ratio e = 0.1, and the scaled mass of the central object is GM, = 1. The two-dimensional pressure and 
plasma beta profile are shown in Fig. [SJ This profile shows that the pressure reaches its maximum within the disk, instead 
of at the inner boundary of the disk, as in the case of a constant density. This difference can be explained as follows. 
Here, also the term po&(R, Z)/T p = p&(R, Z) dominates in the pressure equation (fTO)) . but, due to the small inverse 
aspect ratio, the variation in the gravitational potential is small compared to the one in the density. This implies that the 
density, which is a flux function, mainly determines the shape of the dominant term in the Eq. (fTOj). The temperature 
(not shown) decreases monotonically outwards. The range of plasma beta = [0.082,0.158] and its maximum value is 
reached within the accretion disk. These low values for the plasma beta imply that the disk is strongly magnetized. 
Furthermore, the range of the ratio i/^/vKepier = [0.085, 1.044]. 

For this case, the continuous spectrum has been computed for toroidal mode number n = 1 and poloidal mode 
numbers m — [—3,5]. The results are shown in Fig. which contains two plots. In the left one, the real part of 
the eigenvalues is plotted against the "radial" coordinate s = \ftp, while the right one shows the eigenvalues in the 
complex plane. For this case, there are neither purely exponential unstable or overstable modes. Due to the presence 
of toroidal flow, gravity and a non-constant density, it is possible to create new or wider gaps between the different 
MHD continuum branches, with avoided crossings as a result of poloidal mode couplings. This has been shown by van 
der Hoist et al. (2000b| for tokamak equilibria with toroidal flow. In these gaps, there may exist global modes similar 
to the Toroidal Flow-induced Alfven Eigenmode (TFAE) for tokamaks found by van der Hoist et al. (|2000 a) or to the 
Stratification-induced Alfven Eigenmode (SAE) for coronal flux tubes found by Belien et al. (1996). We will investigate 
this possibility for novel global modes in accretion disk context in future work. 

6.3. Thick accretion disk 

Next, we increase the inverse aspect ratio e from 0.1 to 0.5, while keeping all other quantities the same. Again, the 
poloidal cross-section is chosen circular. The resulting pressure profile is presented in Fig.[8l Here, The pressure decreases 
monotonically outward, similar to the equilibrium with constant density. Again, the term p^>(R,Z) dominates in the 
pressure equation (|13p . The inverse aspect ratio is now not small so that both the gravitational potential and the density 
determines the shape of the pressure profile. Also in this case, the temperature decreases in the outward direction. The 
plasma beta (3 of this accretion disk ranges from 0.237 up to 0.970 and reaches its maximum within the disk. The 
deviation from a Keplerian disk is expressed by the ratio w^/wKcpicr = [0.007,0.282]. Thus the disk completely rotates 
sub-Keplerian. 

For this equilibrium, the continuous MHD spectrum is computed, again for toroidal mode number n = — 1 and poloidal 
mode numbers m — [—3, 5]. A part of the spectra is plotted in Fig. [5] Plot (a) shows the real part of the continuous MHD 
spectrum, while plot (c) shows the imaginary part as a function of the radial flux coordinate s = y/ljj. Plot (c) contains 
the eigenfrequencies plotted in the complex plane. This shows the existence of overstable modes in a thick accretion disk 
due to the presence of toroidal flow and gravity. If one carefully examines the plot (a) one sees that one of the continuum 
branches splits in two at s w 0.87, where the imaginary part becomes zero and at s ~ 0.91 that these two continuum 
branches merge again, whereas the imaginary part becomes non-zero. 

6.4. Maximum growth rate of the MHD continua of thin and thick accretion disks 

The last two cases only differ by the value of the inverse aspect ratio e. We have systematically investigated the influence 
of this ratio on the existence of overstable modes by varying its value, while keeping all other equilibria quantities 
constant. In this parametric study, we essentially move the disk torus radially inwards while keeping the radius of the 
last closed flux surface a constant. We then gradually evolve from the equilibrium shown in Fig. [6] (e = 0.1) to the one 
shown in Fig. [8] (e = 0.5). The growth rate of the most overstable continuum mode is plotted as a function of the inverse 
aspect ratio e in Fig. [TO] Each point of the figure represents an equilibrium computed by FINESSE, from which the MHD 
continua is computed by PHOENIX. In this figure, one can clearly see that no overstable modes are present in these 
disks for e = [0.10, 0.24]. From e = 0.24 up to 0.32, there is a strong increase of the growth rate of the most overstable 
mode, while from e = 0.32 up to at least 0.50 the growth rate remains approximately the same. 
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7. Conclusions 

We have presented the equations describing an axisymmetric MHD accretion disk equilibrium and used FINESSE to find 
numerical solutions. We have considered thin as well as thick disks with a plasma beta of the order 0.1 where the largest 
part of the disk rotates sub-Keplerian. The numerical solutions show that the pressure profile can only decrease in the 
outward direction if the gravitational potential dominates in the pressure equation (|13[) . 

The stability of the disk equilibria have been analysed theoretically as well as numerically. From theory, a stability 
criterion has been derived which closely resembles the Schwarzschild criterion. It is shown that in a disk where the density 
is a flux function, the MHD continua can turn overstable due to the presence of toroidal rotation and gravity. This 
instability is not possible for disks where the temperature or the entropy is a flux function. The numerical results support 
our theoretical conclusion. The effect of the gravity on the overstable modes can either be stabilizing or destabilizing. 
In the case of non-constant density the inverse aspect ratio e has been varied. This study shows that the continua turn 
overstable if the inverse aspect ratio is larger than a certain critical value, approximately 0.24. These localised instabilities 
will lead to small-scale, turbulent dynamics when the equilibria are used in initial conditions for non-linear simulations. 
This dynamics is unrelated the usual magneto-rotational instability. Indeed, these instabilities, as well as the poloidal 
flow driven unstable continuum modes investigated by Goedbloed et al. I2004al provide a new route to MHD turbulence 
in accretion disks. Furthermore, due to the presence of gravity, it is possible to widen or create new gaps in the MHD 
continua. It is possible to show this theoretically by expanding the eigenfunctions in poloidal harmonics. In these gaps 
new types of global ideal MHD eigenmodes can exist. The investigation of these global eigenmodes is left for future work. 
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Fig. 3. (a) Real, and (b) imaginary parts of the sub-spectrum of the MHD continua as functions of the radial flux 
coordinate s = \f^i and (c) represented in the complex plane for toroidal mode number n = — 1 and poloidal mode 
numbers m = [—3, 5]. The corresponding disk equilibrium shown in Fig. [5] has a constant density distribution. 



14 



J. W. S. Blokland et al.: Unstable MHD continuous spectrum of accretion disks 



0.25 




Fig. 4. The growth rate of the most unstable continuum mode as a function of the value of gravitational potential on 
the magnetic axis Rm- The gravitational potential is scaled with respect to the Alfven speed on the magnetic axis. 




Fig. 5. The Brunt- Vaisala frequency projected on a flux surface for the points (a), (b), and (c) shown in Fig. 0] shown 
as a contour plot in a poloidal cross-section. 




Fig. 6. The pressure (gray-scale) and plasma beta (3 = 2p/B 2 (contours) distribution in a poloidal cross-section for an 
accretion disk with density p oc (1 — OAtp) and an inverse aspect ratio e = 0.1. 
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Fig. 7. For the accretion disk equilibrium shown in Fig. (a) the real parts of the sub-spectrum of MHD continua as 
function of the radial flux coordinate s = ^/lp and (b) sub-spectrum of the MHD continua in the complex plane are 
shown for toroidal mode number n — — 1 and poloidal mode numbers m = [—3, 5]. For these mode numbers, no unstable 
or overstable modes are found. 




Fig. 8. The pressure (gray-scale) and plasma beta (3 = 2p/B 2 (contours) for a thick accretion disk with density p oc 
(1 — OAip) and an inverse aspect ratio e = 0.5. 
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Fig. 9. For the thick disk equilibrium shown in Fig. [8l (a) the real and (b) imaginary parts of the sub-spectrum of the 
MHD continua as a function of the radial flux coordinate s = \fip and (c) sub-spectrum of the MHD continua in the 
complex plane are shown for toroidal mode number n = — 1 and poloidal mode numbers m — [—3, 5]. 
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Fig. 10. Growth rate of the most overstable continuum mode as a function of the inverse aspect ratio 



